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ANNE GELB AND EITAN TADMOR 

Abstract. We are concerned with the detection of edges - the location and amplitudes of 
jump discontinuities of piecewise smooth data realized in terms of its discrete grid values. We 
discuss the interplay between two approaches. One approach, realized in the physical space, 
is based on local differences and is typically limited to low-order of accuracy. An alternative 
approach developed in our previous work [6] and realized in the dual Fourier space, is based 
on concentration factors; with a proper choice of concentration factors one can achieve higher- 
orders - in fact in [7] we constructed exponentially accurate edge detectors. Since the stencil 
of these highly-accurate detectors is global, an outside threshold parameter is required to avoid 
oscillations in the immediate neighborhood of discontinuities. In this paper we introduce an 
adaptive edge detection procedure based on a cross-breading between the local and global 
detectors. This is achieved by using the minmod limiter to suppress spurious oscillations 
near discontinuities while retaining high-order accuracy away from the jumps. The resulting 
method provides a family of robust, parameter-free edge-detectors for piecewise smooth data. 
We conclude with a series of one- and two-dimensional simulations. 
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1. Introduction and motivations 

The detection of edges in images is an increasingly important area of research. Not only 
is the knowledge of the jump locations essential in high resolution reconstruction methods, 
e.g. [9, 15, 16], but there are many scientific applications in which the knowledge of the edges 
and their associated jump values are useful in of themselves. One classical example arises in 
magnetic resonance imaging (MRI) segmentation which enables various features of interest to 
be “separated” out for closer examination. Edge detection determines the boundaries of each 
particular region and classifies the type of region based on the jump value at the boundaries. 
Since the most relevant information is often found near the borders of each segmented region, 
it is imperative that edge detection is performed successfully. Furthermore, due to the cost of 
MRI, resolution is somewhat limited. Hence features in an image might be spread over as few as 
two pixel data points with the edges located extremely close together. Finally, computational 
efficiency and robustness are particularly relevant in the case of MRI since huge amounts of 
data must be processed in short periods of time. 

In our previous work, [6, 7, 8], we describe an edge detection procedure which recovers the 
location and amplitudes of edges from spectral information provided either by the continuous 
spectral coefficients or the discrete ones based on equally spaced grid values in physical space. 
These detectors effectively resolve the jump function, [f](x), so that edges are detected by 
separation. Specifically, the detectors “concentrate” near the 0(1) scale of jump discontinuities 
which are effectively separated from the smooth regions where [f](x) ~ 0. This procedure offers 
a large family of edge detectors, each detector is associated with its own particular concentration 
factors. The prescribed procedure is computationally efficient and robust; there are, whoever, 
two major drawbacks. First, in order to “pinpoint” the edges one has to introduce an outside 
threshold parameter to quantify the “large” jumps, i.e., those with [f](x) O(Ax), implying 
that the edge detection method is inherently problem dependent [7]. Second, oscillations form 
in the neighborhood of the jump discontinuities. The particular behavior of these oscillations 
depends on the specific concentration factors used. Therefore it can be difficult to distinguish 
what constitutes a true jump discontinuity as opposed to an oscillating artifact, particularly 
when several jump discontinuities are located in the same neighborhood, i.e. when there is 
limited resolution for the problem. 

In this paper we develop an adaptive, parameter-free edge detection procedure based on the 
nonlinear limiting of low- and high-order concentration factors. Here, we tie in the possibility 
of using “local” edge detectors based on equally spaced grid points to “global” concentration 
methods based on pseudo-spectral coefficients. The paper is organized as follows: In §2 we 
discuss local edge detectors based on equally spaced grid points in physical space. In §3 we 
review the concentration method in dual space and describe low- and high-order concentration 
factors. We also establish the consistency of both, the local and global approaches. An adaptive, 
parameter-free edge detection procedure is introduced in §4: here we use nonlinear limiting 
procedure, based on the so-called minmod limiter, to retain the high-order in smooth regions 
while “limiting” the high-order spurious oscillations in the neighborhoods of the jumps by the 
less oscillatory low-order detectors. Finally, §5 contains numerical examples of two dimensional 
applications of our new adaptive method. 

Although we limit our analysis to periodic piecewise smooth functions on [—7r, 7t) , we note 
that our method is easily adapted for the general intervals and that in two dimensions, edge 
detection is performed one dimension at a time. 
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2. Local detection of edges by (undivided) differences 

Consider a periodic piecewise smooth function f[x) on [— -k, tt) for which we wish to identify 
the points of discontinuity. The corresponding jump function can be defined as [f](x) := 
f(x+) — f(x—), where f(x±) are the right and left side limits of the function at x. Suppose 
we are given discrete grid point values of fix), with fj := f(xj ) on equidistant points Xj = 
—n + jAx, j = 0, • • • , 2N, where Ax = 2 n+i - A j um P discontinuity at x — £ is identihed by its 
enclosed grid cell. Xj < f < Xj + 1 , and is characterized by the asymptotic statement 

f [/] (0 + 0{Ax) for j = : £ e [xj, x j+l ] 

(2-1) A f j+ i := f j+l - J) - \ 

{ 0{ Ax) for j ± j v 

From (2.1), it is clear that the function f(x) experiences a jump discontinuity at every grid 
point Xj, and that the determination of what constitutes a jump is based on an asymptotic 
statement which is inherently a user dependent. It is reasonable, for instance, to assign jump 
values to points corresponding to those x^i’s such that |A/) + i | 0(Ax). Therefore the local 

differences (2.1) can be viewed as an edge detection method of first order. 

To demonstrate the application of (2.1), consider the following examples 


(2.2) f a (x) : = 


:^) 5 , *<o, 


^) 5 , £> 0 . 
7T / 7 


sin (x + 1)' 


fb(x) : = 


sin(^) + 1, -f < x < f, 


sin (x — l) 7 , 


x> f. 


Figure 2.1 displays these piecewise smooth functions on equally spaced grid points. 
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Figures 2.2 and 2.3 demonstrate the convergence of the local first order difference method (2.1) 
to the jump function [f](x). 





0.5 



-0.5 





Figure 2.2: Application of (2.1) for IV = 80 to approximate (a) [f a ](x) and (b) [fb\(x). 


iog|[f a ](x)-Af |i1(2 | 



.. 

-3-2-1 0 1 2 3 x 


log|[f B ](x)-Af j .„ 2 | 

Or 



Figure 2.3: The logarithmic error of the difference formula (2.1) for (a) [f a ](x) and (b) [fb](x), 
with N = 40, 80, and 160. 


While the simple construction of differences (2.1) might suffice for determining the jump dis¬ 
continuities for functions like [f a ](x) and [fb](x), it is not difficult to imagine that distinguishing 
a jump discontinuity from a smooth region might become more complicated for functions with 
higher variation or smaller scales. Therefore we seek an edge detection method that yields 
higher order convergence to zero away from the singular support of the function. 

Following (2.1), we consider higher order local difference formulas of order 2p + 1 to obtain 
faster convergence in the smooth regions of /. For example, the first few higher order difference 
formulas yield 

Al / j+ i := fj+i-fj, 

A ' 3 /j+f := ~fj+ 2 + 3/j+i - 3 fj + fj- 1 , 

A ' 5 /j+ 1 := fj+ 3 ~ 5fj+ 2 + 10/, +1 — 10 fj + 5/,_i — fj- 2 , 


(2.4) 
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Table 1. Short table of coefficients for A 2p+1 / ? - + i. 


2 p + 1 \j + l 

j+5 

j+4 

j+3 

j+2 

j+1 

j 

j-1 

j-2 

j-3 

j-4 

1 





1 

-1 





3 




-1 

3 

-3 

1 




5 



1 

-5 

10 

-10 

5 

-1 



7 


-1 

7 

-21 

35 

-35 

21 

-7 

1 


9 

1 

-9 

36 

-84 

126 

-126 

84 

-36 

9 

-1 


and in general we have 

(2-5) ^ +1 /, + . = E(-h +1 (^)w 

i=~p ' ' 

The coefficients for the first few A 2p+1 /- + i are displayed in table 2.1. Away from the neigh¬ 
borhoods of discontinuities of f(x), the sum (2.5) is of order 0(Ax) 2p+1 . We now wish to 
determine its behavior in the neighborhoods of the jump discontinuities, and then use the re¬ 
sults to develop local edge detectors of odd orders (we omit the even order difference formulas 
which require non-symmetric stencils around x J+ i). 

To motivate our derivation of a local edge detection formula, let us examine the application 
of A 5 / J+ i; the results of examples (2.2) are displayed in figure 2.4. 



Figure 2.4: Application of A 5 / J+ i using N = 80 points to (a) [f a \(x) and (b) [/),] (x). 


A simple computation based on (2.4) yields, in agreement with figure 2.4, 


( 2 . 6 ) 


A ' 5 /,+i 


6[/](£)» if f G [Xj,x j+ i\, 

-4[/](£)> if t e fo-i, Xj) U [x j+l , x j+2 \, 
[/] (0, if € e h- 2 , Xj-i] U [x j+ 2 , Xj +3 ], 
0( Ax) 5 , otherwise. 


It is evident from (2.6) that a local edge detection procedure can be developed from the 
general difference formula (2.5) that can identify both the location and the magnitudes of 
the jump discontinuities of piecewise smooth f(x). Specifically, the (undivided) differences are 
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0{ Ax) 2p+1 in smooth regions, outside the immediate neighborhoods of a jumps, whereas within 
the symmetric 2p-cells neighborhood of such jumps located at, say, x = £, we have 


A 2 P+ lf 


■ 1 ) i ®,p[/](0 + 0(Ax), if 


0( Ax) 


where [11] 


% P = (-o' 


2p+l 

k 


i-l) k = 


Hence we consider the local edge detector 


rt" +1 [/]W = «„, p =( 

Qo , P 2 V 


Letting At —> oo, we have 


kr'im 


(-!)—[/]({), if Xj- r < ( < Ij+i+p, 
Qo , P 

0, otherwise, 


implying that as At —> oo, T^ p+1 [f](x) “concentrates” at the jump discontinuities of f(x) Here 
the index / = traces the cell enclosing the jump so that |£ — Xj±i — Ax/2| < Ax/2. We refer 
to (2.10) as the local concentration property and we note the oscillatory behavior of T^ +1 [f](x) 
is increasing with growing p’s. 

Figures 2.5 and 2.6 display the results for applying the fifth order local edge detection method 
(2.9) to determine the jump functions (2.3). It is clear that the neighborhood of each jump 
discontinuity consists of two cells on each side with 

f [/](£) + 0(Ax), if f G [xj,x j+ 1 ], 

T 5 r xl / \ _ J -|[/](0 + 0(Ax), if £ e [Xj-uXj] U [x j+ 1 ,X j+2 ], 

^ I J[/](0 + O(Aa:), if f e [^-_ 2 , ^-_i] U [x j+2 ,x j+3 ], 

I 0(Ax) 5 i otherwise. 


T 8 o[ f J(x) T so[ f b](x) 



Figure 2.5: Application of Tg Q [f](xj) to approximate (a) [/„.] (x) and (b) [fb]{x). 
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Figure 2.6: The logarithmic error of T^[f](x) for (a) [f a \(x) and (b) [f b ](x), with N = 40,80, 
and 160. 


Figure 2.7 compares the logarithmic errors of the local edge detector (2.9) for various orders p. 
We see that as a direct consequence of (2.10), higher order differencing lead to faster convergence 
away from the jump discontinuities, yet near the discontinuities, first order differencing is more 
advantageous, in the sense that no spurious oscillations are produced. 





-18 7 

-20 L 


iog|[fJ(x)-T 8 ^ , [fKx)| 



' 6 1 
-7 t 


. . 

-3-2-1 0 1 2 3 x 


Figure 2.7: The logarithmic error of T^ +1 [f](xj) using p = 0, • • • , 3 for (a) [ f a ](x ) and (b) 
[fb](x). 


3. Edge detection in the dual space: the global approach 

In [6, 7, 8] we introduced a general family of edge detectors, based on the so-called concen¬ 
tration factors which are implemented in the dual Fourier space. The behavior of the edge 
detection procedure was linked to the type of concentration factor employed. The detectors, in 
the generic case, depend on the Fourier interpolant. Hence they lead to global edge detectors in 
the sense that their stencil involve all the discrete data. At the same time, we show below that 
the local edge detectors (2.9) can be viewed as special cases corresponding to specific trigono¬ 
metric concentration factors. Let us first recall how the global concentration edge detection 
method is formulated (see [6] and [8]). Assume we are given discrete grid point data for a peri¬ 
odic piecewise smooth function f(x ) at Xj = , j — 0, • • • ,2 N. The discrete concentration 







ANNE GELB AND EITAN TADMOR 


detector is defined as 


(3.1) 


N 


TnUKx) :=ni } sgn(k)r 


f\k\Ax 


k=—N 


7T 


fke 


ikx 


where f k are the usual pseudo-spectral coefficients, 


(3.2) 


fk = 


2N 

^ ^ /*/ \ —ikxj 


N+ 1 ^ 
1=0 


f( X j 


and r(sfc) are the discrete concentration factors, 


(3.3) 


r(s k ) 


a(s k )sinc( 


7TSk 

~Y 


\k\Ax 

$k ■ 

7T 


sinc(s) 


sin (s)/s. 


Here cr(s k ) are admissible concentration factors at onr disposal; in [7], we have shown that 
the desired concentration property holds provided the following admissibility requirement is 
fulfilled. 


Corollary 3.1. Let cr(-) G C 2 [0,1] be admissible concentration factor so that ds = 1. 

Then concentration property 

[ [/](£)» l f x = i, 

T N[f]( x ) -»■ < as N ^ oo, 

[ 0 otherwise, 

holds for the edge detector (3.1) associated with the discrete concentration factors r(s k ) = 
a(s k )sinc(- §k). 

To gain a better insight on the behavior of such detectors, we rewrite (3.1) as 


(3.4) 


T T N [f](x) = -Ax 


2 N 

£ 

1=0 


f( X j 


N 

k =1 


(kAx\ sin (fc Ax/2) 


7T 


I 


A; Ax/2 


sin k(x 


x« 


and summation by parts then yields 


(3.5) T T N [f](x) = Ax^(/(x i+ i) - f{xj)) cosk ( x - x j+ 1 / 2 )- 


i=o 


fc=i 


The summation on the right is dominated by the discontinuous cell(s) where \f(xj + f) — f{xj)\ ~ 
0(1) while the contributions of the ’smooth’ cells is negligible due to cancelations of oscillations. 
The precise statement of convergence to the jump function reads [8, §3]: 


(3.6) T^[f](x) = [/](£) cos k(x - g j+ i) + 0(Ax\ log Ax|). 

k =1 

Here f J+ i := x^ + i is the midpoint identifying one discontinuous cell. Observe that such edge 
detectors identify all the jumps by concentrating around the support of the jump function at 
finitely many cells. 
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3.1. Polynomial concentration factors: global edge detectors. Several examples of ad¬ 
missible concentration factors were discussed in [6, 7], We begin our discussion with the poly¬ 
nomial concentration factors, 

(3.7) a 2p+1 (s) :={2p+l)s 2p+1 . 

The corresponding r 2p+ 1 factors are r 2p+ i(s) = cr 2p+1 (s)smc(7rs/2), and (3.5) yields 

~ JL /£-Ar\2p+i cos k(x — x,,i) 

(3.8) ^'[/](s) = A*£(/(*, +1 )-/(*,))£(2p + l)(—)- kAx ■ 

j =0 k =1 

For the first order method, p — 0 in (3.8) reads 

2 N 

(3.9) Tn\J](x) = Ax^2(f(x j+ 1 ) - f(xj))D N (x - x j+ i), 

3= 0 

where D^{y) is the usual Dirichlct sum, D^lpy) = (1 + 2 cos ky) / 2 tt . Hence, (3.9) tells 
us that the discrete concentration kernel associated with the first order polynomial factor (3.7) 
amounts to interpolation of the first-order local differences, f(xj+ 1 ) — fixfi), at the intermedi¬ 
ate grid points, x- + i. In a similar fashion, concentration kernels associated with the higher 
order polynomial factors coincide with higher order derivatives of this interpolant. In fact, an 
equivalent formulation for (3.8) reads 

= (-^(sjvTi)" E (2 P + i)(«0 2, V fa = 

+ k=-N 

(3.10) = (-D^+D^f^ [/](-). 

where gk = 2 W +1 S^o(/( 3: j+ 1 ) — f( x j)) e lkXj+ ^ are the discrete Fourier coefficients of the 
interpolant of the differences we met above. Figures 3.1 and 3.2 display the results for applying 
(3.1) with a fifth order polynomial concentration factor (3.7). 



Figure 3.1: Application of Tff [/] (x) with (3.7) to approximate (a) [f 0 ] (x) and (b) If fix). 
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logiigoo-iaywi 


miMW-T&JMI 



Figure 3.2: The logarithmic error of T^[f](x) with (3.10) for (a) [/ a ](x) and (b) [fb](x), with 
N = 40, 80, and 160. 

We note that resulting edge detectors using the polynomial concentration factors are global in 
the sense that they depend on the global interpolant of /, for T^ p+1 [f](x) = (Tjl [f]Y 2p> (x): the 
corresponding stencil involves all th grid values on the 2ti interval. It is interesting to compare 
these results to the local fifth order edge detection method displayed in figures 2.5 and 2.6. The 
global nature of (3.1) is clearly indicated by the oscillations that exist throughout the domain. 


iog|[f a ](x)-T 8 ’ p,1 (x)| i°g|[U(x)-T 8 r'(x)l 



Figure 3.3: The logarithmic error of T^ p+1 \f](x) with polynomial factor (3.7), p = 0,1, 2, 3, for 
( a ) lfa\(x) and (b) [f b ](x). 

Figure 3.3 compares the logarithmic errors for the concentration method (3.1) with poly¬ 
nomial factors of various order. In the neighborhood of the discontinuities, the first order 
polynomial factor yields the best results, due to the strictly local nature of the method. We 
also note that convergence away from the discontinuities is not necessarily improved for higher 
order polynomial concentration factors, which is due to the global dependence of the polyno¬ 
mial based concentration detectors. Consequently, lower order factors may be more useful in 
calculation. 

3.2. Trigonometric factors: back to local differencing. We now turn to our second family 
of trigonometric concentration factors. As shown in (3.10), the polynomial concentration factors 
(3.7) lead to higher order differentiation of the interpolant of the differences of A[f](xj + 1 ). The 
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closely related trigonometric factors we consider below will amount to higher order differences , 
thus realizing in the dual Fourier space the (undivided) differences we discussed earlier in section 
2. Thus, using the trigonometric factors discussed below we recover the local edge detectors in 
their dual space formulations. Indeed, we will arrive at an alternative approach for determining 
the scaling factor g 0 ,p in (2.9) such that (2.10) is satisfied. 

To demonstrate our point, we consider third-order differences, A 3 /-, i. We use the pseudo- 
spectral Fourier representation 

N 

/(V) = E 

k=-N 


where /*. are the pseudo-spectral coefficients (3.2). After summation by parts, the third-order 
difference formula (2.4) reads 


N 


A7 f+ i = E 

k=-N 


j+ 3 ^ 2N + 1 


2 N 

Y fix^e-^ 1 ( - e ikXj+2 + 3e ikXj+1 - 3e ikXj + e ikXj -^ 


1=0 

2 N N 


(3.11) 


2N + 1 
2 


Y /(^) j k{xi ~ xi) ( - e ik2Ax + 3e ikAx - 3 + e 


—ikAx 


1=0 k=—N 

4 2 N N 


kAx 


2N+1 


Y sin ' 3 sin - x i+i)- 


1=0 k= 1 


Notice that the jump discontinuity is associated with a specific grid point value Xj, whereas in 
the similar global formula (3.4), no grid points are assigned. 

The trigonometric identities 

kAx kAx 

2 sm cos kx l+ i = sm kxi+i — sin kxi , —2 sin sm kx l+ 1 = cos kxi + \ — cos kxi , 

yield 

2 3 AA / X JL LAr 

(3.12) A 3 f(x j+ i) = 2N ^ (/fa+i) - /fo)) 5Z sin 2 _ ^ _cos/;; ^i+| 

1 / n A- 1 

Therefore, the third order (p = 1) local edge detector method can be formulated in the dual 
space as 

(3.13) A 3 f(x j+ ,) = 2f£{f](x j+ ,), 

where according to (3.5), Tff is associated with the concentration factors 03 (s) = Cs4s sin 2 ( 7 rs/ 2 ) 
and c 3 = | in agreement with l/g 0 ,i in (2. 8 ).Thus, Tff [/] (x J+ i) are nothing but the local 
differences T :i [f](x ]+ i). 

At this point it is instructive to show the consistency of the local and global approaches 
for derivation of the same edge detectors methods by computing the scaling factors c p with 
the go,p’s derived in (2.9) from the global perspective. Comparing (3.13) and (3.5), we can 
determine the concentration factor cr^kAx/ir) as 

CT 3 (s) = c 3 2 ssm (—). 


(3.14) 
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The admissibility condition in Corollary 3.1 requires the concentration factors to be normalized 
so that (3.4) holds, f Q l a 3 (s)ds/s = 1 which dictates the value of c 3 as 


(3.15) 


e 3 / 2 2 sin 2 C ) V = 1, 

./n ^ 


22 p+l 

} N + 1 


2 N 


N 


kAx > 


XI (/(^i+i) - /(^)) X sin2p cos m 


x j+\ x i+p- 


1=0 


k =1 


yielding c 3 = |, in agreement with c p = X for p = 1 in (2.8). A similar computation for higher 
order factors leads to 

(3.16) = a, 

Thus, the normalization 

(3.17) 
yields 

(3.18) 


c p I 2 2p sin 2p (^-)ds = 1, 
Jo z 


Cp 


v /7rT(2p+ 1) 


2 2 p p( 2 p+l ) 


2 p 

p 


1 

<? 0 ,p’ 


in agreement with the derivation of the focal edge detectors method (2.9). The general (2p + l)- 
order differencing then corresponds to the trigonometric factors 


(3.19) 


7T S 

o- 2 p + i)(s) = c p 2 2p s sin 2p (y). 


log|[fa](x)-TJf a ](x)} 


-SR 


I I I 1-1 


- 3 - 2-10123 


log|[l b](x) - T s j[fJ(x)} 


\ 

|\t = loc 

-4 

w 

u 

\j V V r 


-5 

1 / y 

: u 


V \ 1 

u 

\ 


.. 


- 3 - 2-10123 


Figure 3.4: The logarithmic error of the fifth order (local) trigonometric vs the (global) poly¬ 
nomial factors T 8 T q [/] (xj) applied to (a) [f a ](x) and (b) [fb](x). 

As exhibited in figure 3.4, the trigonometric edge detectors, (3.19), recover the regions of 
smoothness of a piecewise smooth function more accurately than the polynomial factors (3.7). 
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Figure 3.5: The concentration factors in spectral space using N — 80 for (a) the polynomial 
factors (3.7) and (b) the trigonometric factors (3.19). 

Figure 3.5 compares the analogous concentration factors (3.7) and (3.19) for the polynomial 
and trigonometric edge detectors. As is evident from the comparison, both methods have 
increasing concentration factors on [0,1]. Since both factors reduce the impact of the lower 
order modes while sharpening the contribution of the higher order modes, oscillations are 
inherent in the recovery of the jump function. The steeper slope of the concentration factor 
in the polynomial case results in more prevalent oscillations, suggesting that the application of 
a filter with polynomial edge detection method might effectively reduce the oscillations while 
still yielding high order convergence to zero away from the discontinuities. 



Figure 3.6: (a) Error graph for [f a ](x) applying the concentration method with fifth order 
local, polynomial and filtered polynomial concentration factors, (b) Comparison of the filtered 
polynomial concentration factors (3.20 in spectral space shown for k — 0, • • -40. 


For example, a filtered polynomial concentration factor 
(3.20) ^ 2 p+i( s ) = Const ■ a 2p+ i{s)e~ as \ 

where Const is determined by the property (3.1), can help reduce oscillations and achieve 
higher accuracy away from the jump discontinuities. Figure 3.6 examines the application of 
(3.20) in both physical and spectral space. Here we choose filter parameters a = 32 and k = 8. 
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3.3. Exponential concentration factors. We close our list of examples for dual space de¬ 
tectors with a third family of concentration factors introduced in [7], where we introduced the 
exponential concentration factors, 

i f 1 ~ e -1 

(3.21) <j exp (£) = Const ■ e^- 1 ', Const = J exp{—^~ ——)d£. 

They yield exponential convergence to zero away from the discontinuities with reduced oscil¬ 
lations everywhere except in the neighborhood of the jump discontinuities. We conclude by 
noting that the exponential concentration factor can be seen as the limiting case of high order 
p for the polynomial factors, with additional filtering to reduce the oscillations. It is then more 
powerful than the local methods in terms of convergence away from the discontinuities as well 
as localizing the neighborhood around the discontinuity 

Figure 3.7(a) compares the different concentration factors, (3.7), (3.19) and (3.21). The 
error graph of example (2.2(a)) for the concentration method using the exponential factor is 
exhibited in figure 3.7(b). 



Figure 3.7: (a) Comparison of the exponential (7 = 5), polynomial, and trigonometric concen¬ 
tration factors in spectral space shown for k — 0, • • -40. (b) Error graph for [f a ] (x) applying 
the concentration method with (3.21), 7 = 5 and N = 40,80 and 160. 

It is apparent from the discussion and figures above that the exponential factor both reduces 
oscillations and increases convergence away from a small neighborhood of the jump disconti¬ 
nuity. Both the higher order polynomial factors (3.7) in the global concentration method (3.1) 
and the local edge detection method (2.9) improve convergence rates away from the disconti¬ 
nuities, but decrease the accuracy in the neighborhoods of the jump. As is evident in figures 
2.7 and 3.3, the neighborhoods surrounding the jump discontinuities tend to “spread out” as 
the order of the concentration factor (3.7) or (3.19) increases. The results within the neigh¬ 
borhood are slightly worse for the local edge detection method, as the spreading is dictated 
by (2.10) for the neighborhood around the jump discontinuity. Conversely, the neighborhoods 
around the jump discontinuities are not affected by the application of the concentration method 
using the first order polynomial concentration factor (3.8), or equivalently the first order lo¬ 
cal edge detection method (2.9). In this case, the convergence rate in the neighborhood of the 
discontinuities is O(Ax), but steep gradients elsewhere might be falsely identified as jump loca¬ 
tions, since the convergence outside the neighborhood is also only 0(Ax). If the discontinuities 
are located “far enough” apart, meaning that the jump function can be completely resolved, 
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then the exponential concentration factor is perhaps the most powerful of the methods, since 
T^[f](x) = Tjv[/](x) —> 0 rapidly away from the discontinuities, and the oscillations appear 
only within a small neighborhood of the discontinuities. However, as will be seen in §4, inter¬ 
fering oscillations from “close” discontinuities can make it difficult to accurately recover the 
jump function. 


3.4. Enhanced edge detection. While it is evident that the (local and global) edge detection 
methods converge to the singular support of a piecewise smooth function f(x), the oscillations 
in the neighborhoods of the discontinuities and the various orders of convergence away from 
the discontinuities make it necessary to further enhance the results and “pinpoint” the jump 
discontinuity locations exactly. This has been achieved in [7] by separating the vanishing scales 
in the smooth regions from the 0(1) scales in the neighborhoods of the jump discontinuities. 
Specifically, if {£j}fLi denote the locations of the jump discontinuities of f(x), then for admis¬ 
sible concentration factors in (3.1), the separation of scales is enhanced by computing, for some 
q > 1 at our disposal, 

{ V' /2 ([/]Kj)) a . if i = 

E,,n := Af ?/2 (TJ[/](x))» ~ > 

{ 0(N-i / 2 ), if i # 0. 

leading to the enhanced concentration method, 


(3.22) 


£M/]M) = 


j 

if \E q ,N\ 

1 0. 

if 1 E q , N \ 


< J, 


crit • 


Here J cr n is an 0(1) global threshold parameter signifying the minimal amplitude below which 
jump discontinuities are neglected; after all, almost all grid values experience “jumps” and 
the question is their relative size relative w.r.t. the small scale of 1/N. It is important to 
note that since (2.9) and (3.1) actually detect the neighborhoods of the jump discontinuities, 
0(e(N)), rather than the discontinuities themselves, we use the nonlinear enhancement (3.22) in 
a window of 0(e(N)) to pinpoint the exact locations. Specifically, we determine that the exact 
discontinuities are where the largest amplitude \E q , N | > J crit occur in each small neighborhood 
e(N) of the jump discontinuities. Figure 3.8 show the powerful result of (3.22) in pinpointing 
the jump discontinuities. 


E.oHao^PJM) 

1.5 I- 

1 - 

0.5 ^ 

0 .... . . . 
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- 1 - 5 - 

-2 - • 

.. 
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Figure 3.8: Application of (3.22) (E^[f](x)) to approximate (a) [f a ](x) and (b) [f b \(x). 
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This enhancement procedure requires an outside threshold parameter, reflecting the scaling 
of the function. This becomes an impediment for detecting edges in both small scale problems 
as well as problems with steep gradients and high variation. Additionally, pre-determination 
of the neighborhood value e(N) is also necessary for (3.22), leading to the misidentification of 
jump discontinuities when they are located “too” close together. Next we address how to avoid 
the dependence on the outside threshold parameter. 

4. MinMod edge detection 

As is evident from the results in §2, both the global detectors based on high degree poly¬ 
nomials and, in particular, exponential factors, and local edge detectors based on low degree 
differencing, have the underlying feature that if the convergence rate is fast away from the 
neighborhoods of the discontinuities, then there are more oscillations in the neighborhood of 
the jump discontinuities. Hence it can be difficult to determine the true jump location. This 
is further complicated when jump discontinuities are located near to each other, since the os¬ 
cillations occurring in the neighborhoods of each discontinuity interfere with the true jump 
discontinuities. As noted previously, the first order polynomial edge detection (3.8), or equiva¬ 
lently the first order local edge detection (2.9), does not yield oscillations in the neighborhoods 
of the discontinuities, but has slow convergence away from the discontinuities. On the other 
hand, the exponential concentration factor (3.21) produces rapid convergence to zero away from 
the neighborhoods of discontinuities, but suffers from severe oscillations within the neighbor¬ 
hoods. The loss of monotonicity with the increasing order is, of course, the canonical situation 
in many numerical algorithms. Here we introduce an adaptive edge detection procedure that 
realizes the strengths of both methods. We first observe that away from the jumps, one should 
let the exponentially small factors dominate by taking the smallest (in absolute) value between 
the low-order and high-order detectors. As we approach the jump discontinuity, the high-order 
produce spurious oscillations which should be rejected by the low-order detectors: hence, when 
two values disagree in sign, indicating spurious oscillations, the detectors should be set to zero. 
We end up with the so-called minmod limiter which plays a central role in non-oscillatory re¬ 
construction of high-resolution methods for nonlinear conservation laws, (e.g., [10, 14] and the 
references therein), 

(4.1) Tfi ,nm ° d [f](x) - minmod(f^ II [f](x),fJ} r *'lf](x)), p= 1,2,---, 

where the fc-tuple minmod operation takes the form 

{ s ■ mm(|ai|, |a 2 |,... , \a k \) if sgn(ai) = ... = sgn(a k ) := s 
0 , otherwise. 

Here r exp is the discrete analog of the exponential concentration factors (3.21) and igp+i 
represents the polynomial or trigonometric factors of order (2 p + 1). Although the adaptive 
minmod algorithm is effective for any order polynomial, typically first order would be preferable. 
This adaptive algorithm can be extended to include other concentration factors, e.g. 

(4.2) T% m [f](x) = minmod(T e * P lf]{x ), T poly [f](x), f^ l9 [f](x)), 

for any odd order of polynomial and trigonometric factors. For simplicity, we use the minmod 
algorithm as described in (4.1) with the first and third order polynomial factors (3.7). Figure 
4.1 illustrates the improvement in using (4.1). While some residual oscillations remain in figure 
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4.1(a), it is clear that if the minmod is followed by the nonlinear enhancement procedure (3.22) 
the outside scaling parameter becomes less significant in determining the jump locations, as 
shown in figure 4.1(b). 


minmod(T 8 j[fJ(x)) 


E 80 (minmod(T 8 ^fJ(x)) 


0.25 

0 

-0.25 

-0.5 

-0.75 

-1 

-1.25 

-1.5 


0.25 - 


-0.25 - 
-0.5 1 
-0.75 

-1 

-1.25 

-1.5 


-3 


Figure 4.1: (a) The minmod algorithm (4.1) applied to example 2.2(b) and (b) with the non¬ 
linear enhancement using 80 Fourier modes. 

Furthermore, the minmod, algorithm works even when the jump discontinuities are located 
close together. As an example, consider the jump function and associated jump function ex¬ 
hibited in figure 4.2. 
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Figure 4.2: (a) An example of a piecewise smooth function on IV = 80 points and (b) the jump 
function. 

Figure 4.3 compares the results of application of the concentration method (3.1) with various 
concentration factors and the minmod algorithm, ft is evident that the polynomial factor r 2p+ 1 
does not converge to zero fast enough away from the discontinuities, hence the steep gradi¬ 
ents of the function might be misinterpreted as jump discontinuities. On the other hand, the 
concentration method using r exp causes interfering oscillations in the neighborhoods of the dis¬ 
continuities, making it difficult to determine where the true jumps are. The minmod algorithm 
(4.1) ensures the convergence to the jump function without interference of the oscillations. 
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T [f](x) T[f](x) 




Figure 4.3: Edge detection using first order r po i (dotted), r exp (dashed) and the minmod algo¬ 
rithm (solid) with (a) 80 and (b) 160 points. 


At this point, the application of the enhancement procedure introduced in [7] works directly 
to pinpoint the jump discontinuities, without having to determine a neighborhood parameter 
e(N). Furthermore, one can apply the minimization algorithm discussed in [1] and [2] to reduce 
the 0(Ax) error in the approximation of the amplitude of the jump discontinuity. 


5. Numerical applications 

In this section we describe some two dimensional applications using the minmod edge de¬ 
tection algorithm (4.1). In each case, the adaptive edge detection method is performed one 
dimension at a time and is followed by the enhancement procedure (3.22). 

One important application of edge detection is in the segmentation process associated with 
magnetic resonance imaging (MRI). Specifically, analyst are interested in study particular re¬ 
gions of the brain, and edge detection can be effectively employed to extract or “segment” out 
a region of interest. The classical example of the Shepp-Logan phantom image, displayed in 
figure 5.1, is a typical initial test to determine efficacy of imaging techniques [13]. 



Figure 5.1: Contour plot of the Shepp-Logan brain image. 
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Figure 5.2: The concentration method (3.1) applied to the Shepp-Logan phantom on [128 x 128] 
points using (a) the polynomial concentration factor (3.7) and (b) the exponential concentration 
factor (3.21) at the cross section (x, 0). 


Figure 5.2 exhibits the application of the concentration method (3.1) applied using the con¬ 
centration factors (3.7) and (3.21) for a one dimensional cross section of the image. What 
is evident in figure 5.2 is that the different concentration factors yield different convergence 
patterns in the smooth regions of the image. 



Figure 5.3: The minmod procedure (4.1) applied to the Shepp-Logan brain phantom (a) at the 
cross section (x, 0) and (b) with nonlinear enhancement (3.22). 



Figure 5.4: The nonlinear enhancement procedure (3.22) applied to the Shepp-Logan brain 
phantom image. 
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As exhibited in figures 5.3 and 5.4, application of the minmod algorithm (4.1), further 
enhanced by (3.22), enables complete segmentation of the images. 

Let us now consider an example of a real image, a transversal slice of a human brain provided 
by the Gabrieli Lab [3], exhibited in figure 5.5. 



Figure 5.5: Contour plot of the Gabrieli brain image. 


Since the data is real, random noise is prevalent throughout the image. An outside scaling 
parameter is introduced to help distinguish noise from true image data in the Gabrieli brain 
image figures. A discussion of how the effects of noise is reduced on the concentration method 
(3.1) can be found in [2], Figure 5.6 displays the concentration method (3.1) applied to the 
Gabrieli image on a one dimensional cross section using the first order polynomial (3.7) and 
exponential (3.21) concentration factors. As exhibited in figure 5.7(a), the minmod algorithm 
improves the results by “adapting” the concentration method. The nonlinear enhancement 
procedure (3.22), shown in figures 5.7(b) and 5.8(a), completes the clarification of the distinct 
structures of the Gabrieli image. 



Figure 5.6: The concentration method (3.1) applied to the Gabrieli image using (a) the first 
order polynomial concentration factor (3.7) and (b) the exponential concentration factor (3.21) 
at the cross section (x, 0). 
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Figure 5.7: Edge detection of the Gabrieli image by applying the (a) minmod algorithm (4.1) 
(b) minmod and nonlinear enhancement (3.22) procedures at the cross section (x, 0). 

Finally, we see in figure 5.8(b) that the minmod (4.1) and nonlinear enhancement proce¬ 
dures (3.22) enables high resolution reconstruction of the Gabrieli lab MRI via the Gegenbauer 
reconstruction procedure [9]. 





Figure 5.8: Contour plot of the (a) edge detection of the Gabrieli image by applying the minmod 
and nonlinear enhancement procedures one dimension at a time, (b) Gegenbauer reconstruction 
after the edges have been located. 



6. Conclusion 

In this paper we provide a complete description of edge detection from both the physical 
and spectral data view points. As discussed in §3, the local and global edge detection method 
are consistent in their approaches. The minmod procedure (4.1) provides a computationally 
efficient and robust method to produce edges of images without the need of any outside scaling 
parameters. Thus, the nonlinear enhancement procedure (3.22) works efficiently to “pinpoint” 
the edges of an image. The numerical results given throughout this paper and specifically in 
§5 demonstrate the efficacy of our method, particularly in its ability to enable high resolution 
reconstruction and the segmentation of images. 
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